! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
! test for exponentiation operation of quad complex : complex(16) ^ complex(16)

program test
  integer, parameter :: n = 16
  integer, parameter :: m = 2 * n
  real(16), parameter :: q_tol = 5E-33_16
  integer :: i
  complex(16) :: c1, czero
  complex(16) :: rst1 = (1.58879788522658767867845243213578736_16, 33.8712354897555332621148631597531557_16) **&
                        (666.512059852384793252030932094990168_16, 777.112435235231124235243134588945258_16)
  complex(16) :: rst2 = (2.58879788522658767867812344523213578_16, 1.3387123589755235214562114863159753156_16) **&
                        (-15.6912909810929581908241824108598451_16, -17.1240145105015435439281057182309182_16)
  complex(16) :: rst3 = (2.58879788522658767867845246863213578_16, 1.3387165438352354897556114863159753156_16) **&
                        (0.0_16, 0.0_16)
  complex(16) :: rst4 = (2.58879788522658767867845653243213578_16, 1.3387123548975574654662114863159753156_16) **&
                        (1.0_16, 1.0_16)
  complex(16) :: rst5 = (25.8879788526587678678452432138566578_16, 13.387123548975562114863597531161412856_16) **&
                        (2.0_16, 2.0_16)
  complex(16) :: rst6 = (0.0_16, 0.0_16) ** (0.0_16, 0.0_16)
  complex(16) :: rst7 = (0.0_16, 0.0_16) ** (1.0_16, 0.0_16)
  complex(16) :: rst8 = (0.0_16, 0.0_16) ** (2.0_16, 0.0_16)

  complex(16), parameter :: rst9 = (-156.58879788522658767867845243213578_16, 338.41253312531532421342342141431_16) **&
                                   (77.777777777777777777777777777777777_16, 12.153509430320543502530942309482_16)
  complex(16), parameter :: rst10 = (-223.58879788522658767878845243213578_16, 133.12321542153153214321432153124_16) **&
                                    (-17.777777777777777777777777777777777_16, -69.124319098109201701021009919835_16)
  complex(16), parameter :: rst11 = (1.58879788522658767867845243213578_16, 2.3387123548975562114863159753156_16) **&
                                    (0.0_16, 0.0_16)
  complex(16), parameter :: rst12 = (1.58879788522658767867845243213578_16, 2.3387123548975562114863159753156_16) **&
                                    (1.0_16, 1.0_16)
  complex(16), parameter :: rst13 = (1.58879788522658767867845243213578_16, 2.3387123548975562114863159753156_16) **&
                                    (2.0_16, 2.0_16)
  complex(16), parameter :: rst14 = (0.0_16, 0.0_16) ** (0.0_16, 0.0_16)
  complex(16), parameter :: rst15 = (0.0_16, 0.0_16) ** (1.0_16, 0.0_16)
  complex(16), parameter :: rst16 = (0.0_16, 0.0_16) ** (2.0_16, 0.0_16)
  real(16) :: result(2*n), expect(2*n)

  expect(1) = -4.386018490719466939057955319210772E+0505_16
  expect(2) = 9.154323461024201880262136474129078E+0504_16
  expect(3) = 1.421892410109536649471807454083532E-0004_16
  expect(4) = -1.133763234321218618472305329519878E-0004_16
  expect(5) = 1.0_16
  expect(6) = 0.0_16
  expect(7) = 4.31571255998489129159210697515884793E-0002_16
  expect(8) = 1.80786955261392454822393576648835293_16
  expect(9) = 50.4711560733098985709941399969398_16
  expect(10) = 323.107305231437754064134631999531_16
  expect(11) = 1.0_16
  expect(12) = 0.0_16
  expect(13) = 0.0_16
  expect(14) = 0.0_16
  expect(15) = 0.0_16
  expect(16) = 0.0_16
  expect(17) = -2.137209229739169080121194387899624E+0188_16
  expect(18) = 2.701314279888235794010102976824095E+0189_16
  expect(19) = -1.677713398190515453868739708054252E+0035_16
  expect(20) = 5.966124957584669036467521923878582E+0034_16
  expect(21) = 1.0_16
  expect(22) = 0.0_16
  expect(23) = -0.457183629811155189365045191146104049_16
  expect(24) = 0.964586697959439494454334678478251005_16
  expect(25) = -0.721410626512991567804477279502384_16
  expect(26) = -0.881986495681305897045701792157997_16
  expect(27) = 1.0_16
  expect(28) = 0.0_16
  expect(29) = 0.0_16
  expect(30) = 0.0_16
  expect(31) = 0.0_16
  expect(32) = 0.0_16

  result(1) = rst1%re
  result(2) = rst1%im
  result(3) = rst2%re
  result(4) = rst2%im
  result(5) = rst3%re
  result(6) = rst3%im
  result(7) = rst4%re
  result(8) = rst4%im
  result(9) = rst5%re
  result(10) = rst5%im
  result(11) = rst6%re
  result(12) = rst6%im
  result(13) = rst7%re
  result(14) = rst7%im
  result(15) = rst8%re
  result(16) = rst8%im
  result(17) = rst9%re
  result(18) = rst9%im
  result(19) = rst10%re
  result(20) = rst10%im
  result(21) = rst11%re
  result(22) = rst11%im
  result(23) = rst12%re
  result(24) = rst12%im
  result(25) = rst13%re
  result(26) = rst13%im
  result(27) = rst14%re
  result(28) = rst14%im
  result(29) = rst15%re
  result(30) = rst15%im
  result(31) = rst16%re
  result(32) = rst16%im

  do i = 1, m
    if(expect(i) .eq. 0.0_16) then
      if (result(i) .ne. expect(i)) STOP i
    else
      if (abs((result(i) - expect(i)) / expect(i)) .gt. q_tol) STOP i
    endif
  enddo

  print *, 'PASS'
end
